home *** CD-ROM | disk | FTP | other *** search
/ Personal Computer World 2008 February / PCWFEB08.iso / Software / Resources / Developers / XAMPP 1.5.4 / Windows installer / xampp-win32-1.5.4-installer.exe / xampp / php / pear / Math / Matrix.php < prev    next >
Encoding:
PHP Script  |  2005-12-02  |  45.9 KB  |  1,490 lines
  1. <?php
  2. //
  3. // +----------------------------------------------------------------------+
  4. // | PHP version 4.0                                                      |
  5. // +----------------------------------------------------------------------+
  6. // | Copyright (c) 1997-2001 The PHP Group                                |
  7. // +----------------------------------------------------------------------+
  8. // | This source file is subject to version 2.02 of the PHP license,      |
  9. // | that is bundled with this package in the file LICENSE, and is        |
  10. // | available at through the world-wide-web at                           |
  11. // | http://www.php.net/license/2_02.txt.                                 |
  12. // | If you did not receive a copy of the PHP license and are unable to   |
  13. // | obtain it through the world-wide-web, please send a note to          |
  14. // | license@php.net so we can mail you a copy immediately.               |
  15. // +----------------------------------------------------------------------+
  16. // | Authors: Jesus M. Castagnetto <jmcastagnetto@php.net>                |
  17. // +----------------------------------------------------------------------+
  18. // 
  19. // Matrix definition and manipulation package
  20. // 
  21. // $Id: Matrix.php,v 1.6 2003/04/14 11:27:52 jmcastagnetto Exp $
  22. //
  23.  
  24. require_once 'PEAR.php';
  25. require_once 'Math/Vector/Vector.php';
  26.  
  27. /**
  28.  * Defines a matrix object, conceptualized as an array of arrays such that:
  29.  *
  30.  * <pre>
  31.  * [0][0] [0][1] [0][2] ... [0][M]<br>
  32.  * [1][0] [1][1] [1][2] ... [1][M]<br>
  33.  * ...<br>
  34.  * [N][0] [n][1] [n][2] ... [n][M]<br>
  35.  * </pre>
  36.  *
  37.  * i.e. N rows, M colums
  38.  * 
  39.  * Originally this class was part of NumPHP (Numeric PHP package)
  40.  *
  41.  * @author      Jesus M. Castagnetto <jmcastagnetto@php.net>
  42.  * @access      public
  43.  * @version     1.0
  44.  * @package     Math_Matrix
  45.  */
  46. class Math_Matrix {/*{{{*/
  47.  
  48.     // Properties /*{{{*/
  49.     
  50.     /**
  51.      * Contains the array of arrays defining the matrix
  52.      *
  53.      * @access    private
  54.      * @var     array
  55.      * @see     getData()
  56.      */
  57.     var $_data = null;
  58.  
  59.     /**
  60.      * The number of rows in the matrix
  61.      *
  62.      * @access    private
  63.      * @var     integer
  64.      * @see     getSize()
  65.      */
  66.     var $_num_rows = null;
  67.  
  68.     /**
  69.      * The number of columns in the matrix
  70.      *
  71.      * @access    private
  72.      * @var     integer
  73.      * @see     getSize()
  74.      */
  75.     var $_num_cols = null;
  76.  
  77.     /**
  78.      * The smallest value of all matrix cells
  79.      *
  80.      * @access    private
  81.      * @var     float
  82.      * @see     getMin()
  83.      * @see     getMinMax()
  84.      */
  85.     var $_min = null;
  86.  
  87.     /**
  88.      * The biggest value of all matrix cells
  89.      *
  90.      * @access    private
  91.      * @var     float
  92.      * @see     getMax()
  93.      * @see     getMinMax()
  94.      */
  95.     var $_max = null;
  96.  
  97.     /**
  98.      * A flag indicating if the matrix is square
  99.      * i.e. if $this->_num_cols == $this->_num_rows
  100.      * 
  101.      * @access    private
  102.      * @var     boolean
  103.      * @see     isSquare()
  104.      */
  105.     var $_square = false;
  106.  
  107.  
  108.     /**
  109.      * The Euclidean norm for the matrix: sqrt(sum(e[i][j]^2))
  110.      *
  111.      * @access private
  112.      * @var float
  113.      * @see norm()
  114.      */
  115.     var $_norm = null;
  116.  
  117.     /**
  118.      * The matrix determinant
  119.      *
  120.      * @access private
  121.      * @var float
  122.      * @see determinant()
  123.      */
  124.     var $_det = null;
  125.  
  126.     /**
  127.      * Cutoff error used to test for singular or ill-conditioned matrices
  128.      *
  129.      * @access private
  130.      * @var float
  131.      * @see determinant();
  132.      * @see invert()
  133.      */
  134.     var $_epsilon = 1E-18;
  135.  
  136.     /*}}}*/
  137.  
  138.     /**
  139.      * Constructor for the matrix object
  140.      * 
  141.      * @access  public
  142.      * @param   array   $data
  143.      * @return    object    Math_Matrix
  144.      * @see     $_data
  145.      * @see     setData()
  146.      */
  147.     function Math_Matrix ($data=null) {/*{{{*/
  148.         if (!is_null($data))
  149.             $this->setData($data);
  150.     }/*}}}*/
  151.     
  152.     /**
  153.      * Validates the data and initializes the internal variables (except for the determinant).
  154.      *
  155.      * The validation is performed by by checking that
  156.      * each row (first dimension in the array of arrays)
  157.      * contains the same number of colums (e.g. arrays of the
  158.      * same size)
  159.      *
  160.      * @access  public
  161.      * @param   array   $data   array of arrays to create a matrix object
  162.      * @return    mixed    true on success, a PEAR_Error object otherwise
  163.      */
  164.     function setData($data) {/*{{{*/
  165.         $errorObj = PEAR::raiseError('Invalid data, cannot create/modify matrix');
  166.         if (!is_array($data) || !is_array($data[0])) {
  167.             return $errorObj;
  168.         }
  169.         // check that we got a numeric bidimensional array
  170.         // and that all rows are of the same size
  171.         $nc = count($data[0]);
  172.         $nr = count($data);
  173.         $eucnorm = 0;
  174.         for ($i=0; $i < $nr; $i++) {
  175.             if (count($data[$i]) != $nc) {
  176.                 return $errorObj;
  177.             }
  178.             for ($j=0; $j < $nc; $j++) {
  179.                 if (!is_numeric($data[$i][$j])) {
  180.                     return $errorObj;
  181.                 }
  182.                 $data[$i][$j] = (float) $data[$i][$j];
  183.                 $tmp[] = $data[$i][$j];
  184.                 $eucnorm += $data[$i][$j] * $data[$i][$j];
  185.             }
  186.         }
  187.         $this->_num_rows = $nr;
  188.         $this->_num_cols = $nc;
  189.         $this->_square = ($nr == $nc);
  190.         $this->_min = min($tmp);
  191.         $this->_max = max($tmp);
  192.         $this->_norm = sqrt($eucnorm);
  193.         $this->_data = $data;
  194.         $this->_det = null; // lazy initialization ;-)
  195.         return true;
  196.     }/*}}}*/
  197.  
  198.     /**
  199.      * Returns the array of arrays.
  200.      *
  201.      * @access public
  202.      * @return mixed an array of array of numbers on success, a PEAR_Error otherwise
  203.      */
  204.     function getData () {/*{{{*/
  205.         if ($this->isEmpty()) {
  206.             return PEAR::raiseError('Matrix has not been populated');
  207.         } else {
  208.             return $this->_data;
  209.         }
  210.     }/*}}}*/
  211.  
  212.     /**
  213.      * Checks if the matrix has been initialized.
  214.      *
  215.      * @access public
  216.      * @return boolean TRUE on success, FALSE otherwise
  217.      */
  218.     function isEmpty() {/*{{{*/
  219.         return ( empty($this->_data) || is_null($this->_data) );
  220.     }/*}}}*/
  221.  
  222.  
  223.     /**
  224.      * Returns an array with the number of rows and columns in the matrix
  225.      *
  226.      * @access    public
  227.      * @return    mixed    an array of integers on success, a PEAR_Error object otherwise 
  228.      */
  229.     function getSize() {/*{{{*/
  230.         if ($this->isEmpty())
  231.             return PEAR::raiseError('Matrix has not been populated');
  232.         else
  233.             return array($this->_num_rows, $this->_num_cols);
  234.     }/*}}}*/
  235.  
  236.     /**
  237.      * Checks if it is a square matrix (i.e. num rows == num cols)
  238.      *
  239.      * @access public
  240.      * @return boolean TRUE on success, FALSE otherwise
  241.      */
  242.     function isSquare () {/*{{{*/
  243.         if ($this->isEmpty()) {
  244.             return PEAR::raiseError('Matrix has not been populated');
  245.         } else {
  246.             return $this->_square;
  247.         }
  248.     }/*}}}*/
  249.  
  250.     /**
  251.      * Returns the Euclidean norm of the matrix.
  252.      *
  253.      * Euclidean norm = sqrt( sum( e[i][j]^2 ) )
  254.      * 
  255.      * @access public
  256.      * @return mixed a number on success, a PEAR_Error otherwise
  257.      */
  258.     function norm() {/*{{{*/
  259.         if (!is_null($this->_norm)) {
  260.             return $this->_norm;
  261.         } else {
  262.             return PEAR::raiseError('Uninitialized Math_Matrix object');
  263.         }
  264.     }/*}}}*/
  265.  
  266.     /**
  267.      * Returns a new Math_Matrix object with the same data as the current one
  268.      *
  269.      * @access public
  270.      * @return object Math_Matrix
  271.      */
  272.     function &clone() {/*{{{*/
  273.         return new Math_Matrix($this->_data);
  274.     }/*}}}*/
  275.  
  276.  
  277.     /**
  278.      * Sets the value of the element at (row,col)
  279.      *
  280.      * @param integer $row
  281.      * @param integer $col
  282.      * @param numeric $value
  283.      * @access public
  284.      * @return mixed TRUE on success, a PEAR_Error otherwise
  285.      */
  286.     function setElement($row, $col, $value) {/*{{{*/
  287.         if ($this->isEmpty()) {
  288.             return PEAR::raiseError('Matrix has not been populated');
  289.         }
  290.         if ($row >= $this->_num_rows && $col >= $this->_num_cols) {
  291.             return PEAR::raiseError('Incorrect row and column values');
  292.         }
  293.         if (!is_numeric($value)) {
  294.             return PEAR::raiseError('Incorrect value, expecting a number');
  295.         }
  296.         $this->_data[$row][$col] = $value;
  297.         return true;
  298.     }/*}}}*/
  299.  
  300.     /**
  301.      * Returns the value of the element at (row,col)
  302.      *
  303.      * @param integer $row
  304.      * @param integer $col
  305.      * @access public
  306.      * @return mixed a number on success, a PEAR_Error otherwise
  307.      */
  308.     function getElement($row, $col) {/*{{{*/
  309.         if ($this->isEmpty()) {
  310.             return PEAR::raiseError('Matrix has not been populated');
  311.         }
  312.         if ($row >= $this->_num_rows && $col >= $this->_num_cols) {
  313.             return PEAR::raiseError('Incorrect row and column values');
  314.         }
  315.         return $this->_data[$row][$col];
  316.     }/*}}}*/
  317.  
  318.     /**
  319.      * Returns the row with the given index
  320.      *
  321.      * This method checks that matrix has been initialized and that the
  322.      * row requested is not outside the range of rows.
  323.      *
  324.      * @param integer $row
  325.      * @param optional boolean $asVector whether to return a Math_Vector or a simple array. Default = false.
  326.      * @access public
  327.      * @return mixed an array number on success, a PEAR_Error otherwise
  328.      */
  329.     function getRow ($row, $asVector = false) {/*{{{*/
  330.         if ($this->isEmpty()) {
  331.             return PEAR::raiseError('Matrix has not been populated');
  332.         }
  333.         if (is_integer($row) && $row >= $this->_num_rows) {
  334.             return PEAR::raiseError('Incorrect row value');
  335.         }
  336.         if ($asVector) {
  337.             $classes = get_declared_classes();
  338.             if (!in_array("math_vector", $classes) || !in_array("math_vectopop", $classes)) {
  339.                 return PEAR::raiseError ("Classes Math_Vector and Math_VectorOp undefined". 
  340.                                     " add \"require_once 'Math/Vector/Vector.php'\" to your script");
  341.             }
  342.             return new Math_Vector($this->_data[$row]);
  343.         } else {
  344.             return $this->_data[$row];
  345.         }
  346.     }/*}}}*/
  347.  
  348.     /**
  349.      * Sets the row with the given index to the array
  350.      *
  351.      * This method checks that the row is less than the size of the matrix
  352.      * rows, and that the array size equals the number of columns in the
  353.      * matrix.
  354.      *
  355.      * @param integer $row index of the row
  356.      * @param array $arr array of numbers
  357.      * @access public
  358.      * @return mixed TRUE on success, a PEAR_Error otherwise
  359.      */
  360.     function setRow ($row, $arr) {/*{{{*/
  361.         if ($this->isEmpty()) {
  362.             return PEAR::raiseError('Matrix has not been populated');
  363.         }
  364.         if ($row >= $this->_num_rows) {
  365.             return PEAR::raiseError('Row index out of bounds');
  366.         }
  367.         if (count($arr) != $this->_num_cols) {
  368.             return PEAR::raiseError('Incorrect size for matrix row: expecting '.$this->_num_cols
  369.                     .' columns, got '.count($arr).' columns');
  370.         }
  371.         for ($i=0; $i < $this->_num_cols; $i++) {
  372.             if (!is_numeric($arr[$i])) {
  373.                 return PEAR::raiseError('Incorrect values, expecting numbers');
  374.             }
  375.         }
  376.         $this->_data[$row] = $arr;
  377.         return true;
  378.     }/*}}}*/
  379.  
  380.     /**
  381.      * Returns the column with the given index
  382.      *
  383.      * This method checks that matrix has been initialized and that the
  384.      * column requested is not outside the range of column.
  385.      *
  386.      * @param integer $col
  387.      * @param optional boolean $asVector whether to return a Math_Vector or a simple array. Default = false.
  388.      * @access public
  389.      * @return mixed an array number on success, a PEAR_Error otherwise
  390.      */
  391.     function getCol ($col, $asVector=false) {/*{{{*/
  392.         if ($this->isEmpty()) {
  393.             return PEAR::raiseError('Matrix has not been populated');
  394.         }
  395.         if (is_integer($col) && $col >= $this->_num_cols) {
  396.             return PEAR::raiseError('Incorrect column value');
  397.         }
  398.         for ($i=0; $i < $this->_num_rows; $i++) {
  399.             $ret[$i] = $this->getElement($i,$col);
  400.         }
  401.         if ($asVector) {
  402.             $classes = get_declared_classes();
  403.             if (!in_array("math_vector", $classes) || !in_array("math_vectopop", $classes)) {
  404.                 return PEAR::raiseError ("Classes Math_Vector and Math_VectorOp undefined". 
  405.                                     " add \"require_once 'Math/Vector/Vector.php'\" to your script");
  406.             }
  407.             return new Math_Vector($ret);
  408.         } else {
  409.             return $ret;
  410.         }
  411.     }/*}}}*/
  412.  
  413.     /**
  414.      * Sets the column with the given index to the array
  415.      *
  416.      * This method checks that the column is less than the size of the matrix
  417.      * columns, and that the array size equals the number of rows in the
  418.      * matrix.
  419.      *
  420.      * @param integer $col index of the column
  421.      * @param array $arr array of numbers
  422.      * @access public
  423.      * @return mixed TRUE on success, a PEAR_Error otherwise
  424.      */
  425.     function setCol ($col, $arr) {/*{{{*/
  426.         if ($this->isEmpty()) {
  427.             return PEAR::raiseError('Matrix has not been populated');
  428.         }
  429.         if ($col >= $this->_num_cols) {
  430.             return PEAR::raiseError('Incorrect column value');
  431.         }
  432.         if (count($arr) != $this->_num_cols) {
  433.             return PEAR::raiseError('Incorrect size for matrix column');
  434.         }
  435.         for ($i=0; $i < $this->_num_rows; $i++) {
  436.             if (!is_numeric($arr[$i])) {
  437.                 return PEAR::raiseError('Incorrect values, expecting numbers');
  438.             } else {
  439.                 $err = $this->setElement($i, $col, $arr[$i]);
  440.                 if (PEAR::isError($err)) {
  441.                     return $err;
  442.                 }
  443.             }
  444.             
  445.         }
  446.         return true;
  447.     }/*}}}*/
  448.  
  449.     /**
  450.      * Swaps the rows with the given indices
  451.      *
  452.      * @param integer $i
  453.      * @param integer $j
  454.      * @access public
  455.      * @return mixed TRUE on success, a PEAR_Error otherwise
  456.      */
  457.     function swapRows($i, $j) {/*{{{*/
  458.         $r1 = $this->getRow($i);
  459.         if (PEAR::isError($r1)) {
  460.             return $r1;
  461.         }
  462.         $r2 = $this->getRow($j);
  463.         if (PEAR::isError($r2)) {
  464.             return $r2;
  465.         }
  466.         $e = $this->setRow($j, $r1);
  467.         if (PEAR::isError($e)) {
  468.             return $e;
  469.         }
  470.         $e = $this->setRow($i, $r2);
  471.         if (PEAR::isError($e)) {
  472.             return $e;
  473.         }
  474.         return true;
  475.     }/*}}}*/
  476.  
  477.     /**
  478.      * Swaps the columns with the given indices
  479.      *
  480.      * @param integer $i
  481.      * @param integer $j
  482.      * @access public
  483.      * @return mixed TRUE on success, a PEAR_Error otherwise
  484.      */
  485.     function swapCols($i, $j) {/*{{{*/
  486.         $r1 = $this->getCol($i);
  487.         if (PEAR::isError($r1)) {
  488.             return $r1;
  489.         }
  490.         $r2 = $this->getCol($j);
  491.         if (PEAR::isError($r2)) {
  492.             return $r2;
  493.         }
  494.         $e = $this->setCol($j, $r1);
  495.         if (PEAR::isError($e)) {
  496.             return $e;
  497.         }
  498.         $e = $this->setCol($i, $r2);
  499.         if (PEAR::isError($e)) {
  500.             return $e;
  501.         }
  502.         return true;
  503.     }/*}}}*/
  504.  
  505.     /**
  506.      * Swaps a given row with a given column. Only valid for square matrices.
  507.      *
  508.      * @param integer $row index of row
  509.      * @param integer $col index of column
  510.      * @access public
  511.      * @return mixed TRUE on success, a PEAR_Error otherwise
  512.      */
  513.     function swapRowCol ($row, $col) {/*{{{*/
  514.         if (!$this->isSquare() || !is_int($row) || !is_int($col)) {
  515.             return PEAR::raiseError("Parameters must be row and a column indices");
  516.         }
  517.         $c = $this->getCol($col);
  518.         if (PEAR::isError($c)) {
  519.             return $c;
  520.         }
  521.         $r = $this->getRow($row);
  522.         if (PEAR::isError($r)) {
  523.             return $r;
  524.         }
  525.         $e = $this->setCol($col, $r);
  526.         if (PEAR::isError($e)) {
  527.             return $e;
  528.         }
  529.         $e = $this->setRow($row, $c);
  530.         if (PEAR::isError($e)) {
  531.             return $e;
  532.         }
  533.         return true;
  534.     }/*}}}*/
  535.  
  536.     /**
  537.      * Returns the minimum value of the elements in the matrix
  538.      *
  539.      * @access public
  540.      * @return mixed a number on success, a PEAR_Error otherwise
  541.      */
  542.     function getMin () {/*{{{*/
  543.         if ($this->isEmpty()) {
  544.             return PEAR::raiseError('Matrix has not been populated');
  545.         } else {
  546.             return $this->_min;
  547.         }
  548.     }/*}}}*/
  549.     
  550.     /**
  551.      * Returns the maximum value of the elements in the matrix
  552.      *
  553.      * @access public
  554.      * @return mixed a number on success, a PEAR_Error otherwise
  555.      */
  556.     function getMax () {/*{{{*/
  557.         if ($this->isEmpty()) {
  558.             return PEAR::raiseError('Matrix has not been populated');
  559.         } else {
  560.             return $this->_max;
  561.         }
  562.     }/*}}}*/
  563.  
  564.     /**
  565.      * Gets the position of the first element with the given value
  566.      *
  567.      * @param numeric $val
  568.      * @access public
  569.      * @return mixed an array of two numbers on success, FALSE if value is not found, and PEAR_Error otherwise
  570.      */
  571.     function getValueIndex ($val) {/*{{{*/
  572.         if ($this->isEmpty()) {
  573.             return PEAR::raiseError('Matrix has not been populated');
  574.         }
  575.         for ($i=0; $i < $this->_num_rows; $i++) {
  576.             for ($j=0; $j < $this->_num_cols; $j++) {
  577.                 if ($this->_data[$i][$j] == $val) {
  578.                     return array($i, $j);
  579.                 }
  580.             }
  581.         }
  582.         return false;
  583.     }/*}}}*/
  584.  
  585.     /**
  586.      * Gets the position of the element with the minimum value
  587.      *
  588.      * @access public
  589.      * @return mixed an array of two numbers on success, FALSE if value is not found, and PEAR_Error otherwise
  590.      * @see getValueIndex()
  591.      */
  592.     function getMinIndex () {/*{{{*/
  593.         if ($this->isEmpty()) {
  594.             return PEAR::raiseError('Matrix has not been populated');
  595.         } else {
  596.             return $this->getValueIndex($this->_min);
  597.         }
  598.     }/*}}}*/
  599.  
  600.     /**
  601.      * Gets the position of the element with the maximum value
  602.      *
  603.      * @access public
  604.      * @return mixed an array of two numbers on success, FALSE if value is not found, and PEAR_Error otherwise
  605.      * @see getValueIndex()
  606.      */
  607.     function getMaxIndex () {/*{{{*/
  608.         if ($this->isEmpty()) {
  609.             return PEAR::raiseError('Matrix has not been populated');
  610.         } else {
  611.             return $this->getValueIndex($this->_max);
  612.         }
  613.     }/*}}}*/
  614.  
  615.     /**
  616.      * Transpose the matrix rows and columns
  617.      *
  618.      * @access public
  619.      * @return mixed TRUE on success, PEAR_Error otherwise
  620.      */
  621.     function transpose () {/*{{{*/
  622.         if (!$this->isSquare()) {
  623.             return PEAR::raiseError("Transpose is undefined for non-sqaure matrices");
  624.         }
  625.         list($nr, $nc) = $this->getSize();
  626.         $data = array();
  627.         for ($i=0; $i < $nc; $i++) {
  628.             $col = $this->getCol($i);
  629.             if (PEAR::isError($col)) {
  630.                 return $col;
  631.             } else {
  632.                 $data[] = $col;
  633.             }
  634.         }
  635.         return $this->setData($data);
  636.     }/*}}}*/
  637.  
  638.     /**
  639.      * Returns the trace of the matrix. Trace = sum(e[i][j]), for all i == j
  640.      *
  641.      * @access public
  642.      * @return mixed a number on success, PEAR_Error otherwise
  643.      */
  644.     function trace() {/*{{{*/
  645.         if ($this->isEmpty()) {
  646.             return PEAR::raiseError('Matrix has not been populated');
  647.         }
  648.         if (!$this->isSquare()) {
  649.             return PEAR::raiseError('Trace undefined for non-square matrices');
  650.         }
  651.         $trace = 0;
  652.         for ($i=0; $i < $this->_num_rows; $i++) {
  653.             $trace += $this->getElement($i, $i);
  654.         }
  655.         return $trace;
  656.     }/*}}}*/
  657.     
  658.     /**
  659.      * Calculates the matrix determinant using Gaussian elimination with partial pivoting.
  660.      *
  661.      * At each step of the pivoting proccess, it checks that the normalized
  662.      * determinant calculated so far is less than 10^-18, trying to detect
  663.      * singular or ill-conditioned matrices
  664.      *
  665.      * @access public
  666.      * @return mixed a number on success, a PEAR_Error otherwise
  667.      */
  668.     function determinant() {/*{{{*/
  669.         if (!is_null($this->_det) && is_numeric($this->_det)) {
  670.             return $this->_det;
  671.         }
  672.         if ($this->isEmpty()) {
  673.             return PEAR::raiseError('Matrix has not been populated');
  674.         }
  675.         if (!$this->isSquare()) {
  676.             return PEAR::raiseError('Determinant undefined for non-square matrices');
  677.         }
  678.         $norm = $this->norm();
  679.         if (PEAR::isError($norm)) {
  680.             return $norm;
  681.         }
  682.         $det = 1.0;
  683.         $sign = 1;
  684.         // work on a copy
  685.         $m = $this->clone();
  686.         list($nr, $nc) = $m->getSize();
  687.         for ($r=0; $r<$nr; $r++) { 
  688.             // find the maximum element in the column under the current diagonal element
  689.             $ridx = $m->_maxElementIndex($r);
  690.             if (PEAR::isError($ridx)) {
  691.                 return $ridx;
  692.             }
  693.             if ($ridx != $r) {
  694.                 $sign = -$sign;
  695.                 $e = $m->swapRows($r, $ridx);
  696.                 if (PEAR::isError($e)) {
  697.                     return $e;
  698.                 }
  699.             }
  700.             // pivoting element
  701.             $pelement = $m->getElement($r, $r);
  702.             if (PEAR::isError($pelement)) {
  703.                 return $pelement;
  704.             }
  705.             $det *= $pelement;
  706.             // Is this an singular or ill-conditioned matrix?
  707.             // i.e. is the normalized determinant << 1 and -> 0?
  708.             if ((abs($det)/$norm) < $this->_epsilon) {
  709.                 return PEAR::raiseError('Probable singular or ill-conditioned matrix, normalized determinant = '
  710.                         .(abs($det)/$norm));
  711.             }
  712.             if ($pelement == 0) {
  713.                 return PEAR::raiseError('Cannot continue, pivoting element is zero');
  714.             }
  715.             // zero all elements in column below the pivoting element
  716.             for ($i = $r + 1; $i < $nr; $i++) {
  717.                 $factor = $m->getElement($i, $r) / $pelement;
  718.                 for ($j=$r; $j < $nc; $j++) {
  719.                     $val = $m->getElement($i, $j) - $factor*$m->getElement($r, $j);
  720.                     $e = $m->setElement($i, $j, $val);
  721.                     if (PEAR::isError($e)) {
  722.                         return $e;
  723.                     }
  724.                 }
  725.             }
  726.             // for debugging
  727.             //echo "COLUMN: $r\n";
  728.             //echo $m->toString()."\n";
  729.         }
  730.         unset($m);
  731.         if ($sign < 0) {
  732.             $det = -$det;
  733.         }
  734.         // save the value
  735.         $this->_det = $det;
  736.         return $det;
  737.     }/*}}}*/
  738.  
  739.     /**
  740.      * Returns the normalized determinant = abs(determinant)/(euclidean norm)
  741.      *
  742.      * @access public
  743.      * @return mixed a positive number on success, a PEAR_Error otherwise
  744.      */
  745.     function normalizedDeterminant() {/*{{{*/
  746.         $det = $this->determinant();
  747.         if (PEAR::isError($det)) {
  748.             return $det;
  749.         }
  750.         $norm = $this->norm();
  751.         if (PEAR::isError($norm)) {
  752.             return $norm;
  753.         }
  754.         if ($norm == 0) {
  755.             return PEAR::raiseError('Undefined normalized determinant, euclidean norm is zero');
  756.         }
  757.         return abs($det / $norm);
  758.     }/*}}}*/
  759.  
  760.     /**
  761.      * Returns the index of the row with the maximum value under column of the element e[i][i]
  762.      *
  763.      * @access protected
  764.      * @return an integer
  765.      */
  766.     function _maxElementIndex($r) {/*{{{*/
  767.         $max = 0;
  768.         $idx = -1;
  769.         list($nr, $nc) = $this->getSize();
  770.         $arr = array();
  771.         for ($i=$r; $i<$nr; $i++) {
  772.             $val = abs($this->_data[$i][$r]);
  773.             if ($val > $max) {
  774.                 $max = $val;
  775.                 $idx = $i;
  776.             }
  777.         }
  778.         if ($idx == -1) {
  779.             $idx = $r;
  780.         }
  781.         return $idx;
  782.     }/*}}}*/
  783.  
  784.     /**
  785.      * Inverts a matrix using Gauss-Jordan elimination with partial pivoting
  786.      *
  787.      * @access public
  788.      * @return mixed the value of the matrix determinant on success, PEAR_Error otherwise
  789.      * @see scaleRow()
  790.      */
  791.     function invert() {
  792.         if ($this->isEmpty()) {
  793.             return PEAR::raiseError('Matrix has not been populated');
  794.         }
  795.         if (!$this->isSquare()) {
  796.             return PEAR::raiseError('Determinant undefined for non-square matrices');
  797.         }
  798.         $norm = $this->norm();
  799.         $sign = 1;
  800.         $det = 1.0;
  801.         // work on a copy to be safe
  802.         $m = $this->clone();
  803.         list($nr, $nc) = $m->getSize();
  804.         // Unit matrix to use as target
  805.         $q = Math_Matrix::makeUnit($nr);
  806.         for ($i=0; $i<$nr; $i++) {
  807.             $ridx = $this->_maxElementIndex($i);
  808.             if ($i != $ridx) { 
  809.                 $sign = -$sign;
  810.                 $e = $m->swapRows($i, $ridx);
  811.                 if (PEAR::isError($e)) {
  812.                     return $e;
  813.                 }
  814.                 $e = $q->swapRows($i, $ridx);
  815.                 if (PEAR::isError($e)) {
  816.                     return $e;
  817.                 }
  818.             }
  819.             $pelement = $m->getElement($i, $i);
  820.             if (PEAR::isError($pelement)) {
  821.                 return $pelement;
  822.             }
  823.             if ($pelement == 0) {
  824.                 return PEAR::raiseError('Cannot continue inversion, pivoting element is zero');
  825.             }
  826.             $det *= $pelement;
  827.             if ((abs($det)/$norm) < $this->_epsilon) {
  828.                 return PEAR::raiseError('Probable singular or ill-conditioned matrix, normalized determinant = '
  829.                         .(abs($det)/$norm));
  830.             }
  831.             $e = $m->scaleRow($i, 1/$pelement);
  832.             if (PEAR::isError($e)) {
  833.                 return $e;
  834.             }
  835.             $e = $q->scaleRow($i, 1/$pelement);
  836.             if (PEAR::isError($e)) {
  837.                 return $e;
  838.             }
  839.             // zero all column elements execpt for the current one
  840.             $tpelement = $m->getElement($i, $i);
  841.             for ($j=0; $j<$nr; $j++) {
  842.                 if ($j == $i) {
  843.                     continue;
  844.                 }
  845.                 $factor = $m->getElement($j, $i) / $tpelement;
  846.                 for ($k=0; $k<$nc; $k++) {
  847.                     $vm = $m->getElement($j, $k) - $factor * $m->getElement($i, $k);
  848.                     $vq = $q->getElement($j, $k) - $factor * $q->getElement($i, $k);
  849.                     $m->setElement($j, $k, $vm);
  850.                     $q->setElement($j, $k, $vq);
  851.                 }
  852.             }
  853.             // for debugging
  854.             /*
  855.             echo "COLUMN: $i\n";
  856.             echo $m->toString()."\n";
  857.             echo $q->toString()."\n";
  858.             */
  859.         }
  860.         $data = $q->getData();
  861.         /*
  862.         // for debugging
  863.         echo $m->toString()."\n";
  864.         echo $q->toString()."\n";
  865.         */
  866.         unset($m);
  867.         unset($q);
  868.         $e = $this->setData($data);
  869.         if (PEAR::isError($e)) {
  870.             return $e;
  871.         }
  872.         if ($sign < 0) {
  873.             $det = -$det;
  874.         }
  875.         $this->_det = $det;
  876.         return $det;
  877.     }
  878.  
  879.     /**
  880.      * Returns a submatrix from the position (row, col), with nrows and ncols
  881.      *
  882.      * @access public
  883.      * @return object Math_Matrix on success, PEAR_Error otherwise
  884.      */
  885.     function &getSubMatrix ($row, $col, $nrows, $ncols) {/*{{{*/
  886.         if (!is_numeric($row) || !is_numeric($col)
  887.             || !is_numeric($nrows) || !is_numeric($ncols)) {
  888.             return PEAR::raiseError('Parameters must be a initial row and column, and number of rows and columns in submatrix');
  889.         }
  890.         list($nr, $nc) = $this->getSize();
  891.         if ($rows + $nrows > $nr) {
  892.             return PEAR::raiseError('Rows in submatrix more than in original matrix');
  893.         }
  894.         if ($cols + $ncols > $nr) {
  895.             return PEAR::raiseError('Columns in submatrix more than in original matrix');
  896.         }
  897.         $data = array();
  898.         for ($i=0; $i < $nrows; $i++) {
  899.             for ($j=0; $j < $ncols; $j++) {
  900.                 $data[$i][$j] = $this->getElement($i + $row, $j + $col);
  901.             }
  902.         }
  903.         return new Math_Matrix($data);
  904.     }/*}}}*/
  905.     
  906.     /**
  907.      * Returns a simple string representation of the matrix
  908.      *
  909.      * @param optional string $format a sprintf() format used to print the matrix elements (default = '%6.2f')
  910.      * @return mixed a string on success, PEAR_Error otherwise
  911.      */
  912.     function toString ($format='%6.2f') {/*{{{*/
  913.         if ($this->isEmpty()) {
  914.             return PEAR::raiseError('Matrix has not been populated');
  915.         }
  916.         $out = "";
  917.         for ($i=0; $i < $this->_num_rows; $i++) {
  918.             for ($j=0; $j < $this->_num_cols; $j++) {
  919.                 $out .= sprintf($format, $this->_data[$i][$j]);
  920.             }
  921.             $out .= "\n";
  922.         }
  923.         return $out;
  924.     }/*}}}*/
  925.     
  926.     /**
  927.      * Returns an HTML table representation of the matrix elements
  928.      *
  929.      * @access public
  930.      * @return a string on success, PEAR_Error otherwise
  931.      */
  932.     function toHTML() {/*{{{*/
  933.         if ($this->isEmpty()) {
  934.             return PEAR::raiseError('Matrix has not been populated');
  935.         }
  936.         $out = "<table border>\n\t<caption align=\"top\"><b>Matrix</b>";
  937.         $out .= "</caption>\n\t<tr align=\"center\">\n\t\t<th>";
  938.         $out .= $this->_num_rows."x".$this->_num_cols."</th>";
  939.         for ($i=0; $i < $this->_num_cols; $i++) {
  940.             $out .= "<th>".$i."</th>";
  941.         }
  942.         $out .= "\n\t</tr>\n";
  943.         for ($i=0; $i < $this->_num_rows; $i++) {
  944.             $out .= "\t<tr align=\"center\">\n\t\t<th>".$i."</th>";
  945.             for ($j=0; $j < $this->_num_cols; $j++) {
  946.                 $out .= "<td bgcolor=\"#ffffdd\">".$this->_data[$i][$j]."</td>";
  947.             }
  948.             $out .= "\n\t</tr>";
  949.         }
  950.         return $out."\n</table>\n";
  951.     }/*}}}*/
  952.  
  953.  
  954.     // Binary operations
  955.     
  956.     /**
  957.      * Adds a matrix to this one
  958.      *
  959.      * @param object Math_Matrix $m1
  960.      * @return mixed TRUE on success, PEAR_Error otherwise
  961.      * @see getSize()
  962.      * @see getElement()
  963.      * @see setData()
  964.      */
  965.     function add ($m1) {/*{{{*/
  966.         if (!Math_Matrix::isMatrix($m1)) {
  967.             return PEAR::raiseError("Parameter must be a Math_Matrix object");
  968.         }
  969.         if ($this->getSize() != $m1->getSize()) {
  970.             return PEAR::raiseError("Matrices must have the same dimensions");
  971.         }
  972.         list($nr, $nc) = $m1->getSize();
  973.         $data = array();
  974.         for ($i=0; $i < $nr; $i++) {
  975.             for ($j=0; $j < $nc; $j++) {
  976.                 $el1 = $m1->getElement($i,$j);
  977.                 if (PEAR::isError($el1)) {
  978.                     return $el1;
  979.                 }
  980.                 $el = $this->getElement($i,$j);
  981.                 if (PEAR::isError($el)) {
  982.                     return $el;
  983.                 }
  984.                 $data[$i][$j] = $el + $el1;
  985.             }
  986.         }
  987.         if (!empty($data)) {
  988.             return $this->setData($data);
  989.         } else {
  990.             return PEAR::raiseError('Undefined error');
  991.         }
  992.     }/*}}}*/
  993.  
  994.     /**
  995.      * Substracts a matrix from this one
  996.      *
  997.      * @param object Math_Matrix $m1
  998.      * @return mixed TRUE on success, PEAR_Error otherwise
  999.      * @see getSize()
  1000.      * @see getElement()
  1001.      * @see setData()
  1002.      */
  1003.     function sub (&$m1) {/*{{{*/
  1004.         if (!Math_Matrix::isMatrix($m1)) {
  1005.             return PEAR::raiseError("Parameter must be a Math_Matrix object");
  1006.         }
  1007.         if ($this->getSize() != $m1->getSize()) {
  1008.             return PEAR::raiseError("Matrices must have the same dimensions");
  1009.         }
  1010.         list($nr, $nc) = $m1->getSize();
  1011.         $data = array();
  1012.         for ($i=0; $i < $nr; $i++) {
  1013.             for ($j=0; $j < $nc; $j++) {
  1014.                 $el1 = $m1->getElement($i,$j);
  1015.                 if (PEAR::isError($el1)) {
  1016.                     return $el1;
  1017.                 }
  1018.                 $el = $this->getElement($i,$j);
  1019.                 if (PEAR::isError($el)) {
  1020.                     return $el;
  1021.                 }
  1022.                 $data[$i][$j] = $el - $el1;
  1023.             }
  1024.         }
  1025.         if (!empty($data)) {
  1026.             return $this->setData($data);
  1027.         } else {
  1028.             return PEAR::raiseError('Undefined error');
  1029.         }
  1030.     }/*}}}*/
  1031.  
  1032.     /**
  1033.      * Scales the matrix by a given factor
  1034.      *
  1035.      * @param numeric $scale the scaling factor
  1036.      * @return mixed TRUE on success, PEAR_Error otherwise
  1037.      * @see getSize()
  1038.      * @see getElement()
  1039.      * @see setData()
  1040.      */
  1041.     function scale ($scale) {/*{{{*/
  1042.         if (!is_numeric($scale)) {
  1043.             return PEAR::raiseError("Parameter must be a number");
  1044.         }
  1045.         list($nr, $nc) = $this->getSize();
  1046.         $data = array();
  1047.         for ($i=0; $i < $nr; $i++) {
  1048.             for ($j=0; $j < $nc; $j++) {
  1049.                 $data[$i][$j] = $scale * $this->getElement($i,$j);
  1050.             }
  1051.         }
  1052.         if (!empty($data)) {
  1053.             return $this->setData($data);
  1054.         } else {
  1055.             return PEAR::raiseError('Undefined error');
  1056.         }
  1057.     }/*}}}*/
  1058.  
  1059.     /**
  1060.      * Multiplies (scales) a row by the given factor
  1061.      *
  1062.      * @access public
  1063.      * @param integer $row the row index
  1064.      * @param numeric $factor the scaling factor
  1065.      * @return mixed TRUE on success, a PEAR_Error otherwise
  1066.      * @see invert()
  1067.      */
  1068.     function scaleRow($row, $factor) {/*{{{*/
  1069.         if ($this->isEmpty()) {
  1070.             return PEAR::raiseError('Uninitialized Math_Matrix object');
  1071.         }
  1072.         if (!is_integer($row) || !is_numeric($factor)) {
  1073.             return PEAR::raiseError('Row index must be an integer, and factor a valid number');
  1074.         }
  1075.         if ($row >= $this->_num_rows) {
  1076.             return PEAR::raiseError('Row index out of bounds');
  1077.         }
  1078.         $r = $this->getRow($row);
  1079.         if (PEAR::isError($r)) {
  1080.             return $r;
  1081.         }
  1082.         $nr = count($r);
  1083.         for ($i=0; $i<$nr; $i++) {
  1084.             $r[$i] *= $factor;
  1085.         }
  1086.         return $this->setRow($row, $r);
  1087.     }/*}}}*/
  1088.  
  1089.     /**
  1090.      * Multiplies a matrix by this one
  1091.      *
  1092.      * @param object Math_Matrix $m1
  1093.      * @return mixed TRUE on success, PEAR_Error otherwise
  1094.      * @see getSize()
  1095.      * @see getRow()
  1096.      * @see getCol()
  1097.      * @see setData()
  1098.      */
  1099.     function multiply(&$m1) {/*{{{*/
  1100.         $epsilon = 1E-10;
  1101.         if (!Math_Matrix::isMatrix($m1)) {
  1102.             return PEAR::raiseError ('Wrong parameter, expected a Math_Matrix object');
  1103.         }
  1104.         list($nr, $nc) = $this->getSize();
  1105.         list($nr1, $nc1) = $m1->getSize();
  1106.         if ($nc1 != $nr) {
  1107.             return PEAR::raiseError('Incompatible sizes columns in matrix must be the same as rows in parameter matrix');
  1108.         }
  1109.         $data = array();
  1110.         for ($i=0; $i < $nr; $i++) {
  1111.             for ($j=0; $j < $nc1; $j++) {
  1112.                 $data[$i][$j] = 0;
  1113.                 for ($k=0; $k < $nc; $k++) {
  1114.                     $data[$i][$j] += $this->getElement($i,$k) * $m1->getElement($k, $j);
  1115.                 }
  1116.                 // take care of some round-off errors
  1117.                 if ($data[$i][$j] <= $epsilon) {
  1118.                     $data[$i][$j] = 0.0;
  1119.                 }
  1120.             }
  1121.         }
  1122.         if (!empty($data)) {
  1123.             return $this->setData($data);
  1124.         } else {
  1125.             return PEAR::raiseError('Undefined error');
  1126.         }
  1127.     }/*}}}*/
  1128.  
  1129.     /**
  1130.      * Multiplies a vector by this matrix
  1131.      *
  1132.      * @param object Math_Vector $v1
  1133.      * @return object an instance of Math_Vector on success, PEAR_Error otherwise
  1134.      * @see getSize()
  1135.      * @see getRow()
  1136.      * @see Math_Vector::get()
  1137.      */
  1138.     function &vectorMultiply(&$v1) {/*{{{*/
  1139.         if (!Math_VectorOp::isVector($v1)) {
  1140.             return PEAR::raiseError ("Wrong parameter, a Math_Vector object");
  1141.         }
  1142.         list($nr, $nc) = $this->getSize();
  1143.         $nv = $v1->size();
  1144.         if ($nc != $nv) {
  1145.             return PEAR::raiseError("Incompatible number of columns in matrix ($nc) must ".
  1146.                         "be the same as the number of elements ($nv) in the vector");
  1147.         }
  1148.         $data = array();
  1149.         for ($i=0; $i < $nr; $i++) {
  1150.             $data[$i] = 0;
  1151.             for ($j=0; $j < $nv; $j++) {
  1152.                 $data[$i] += $this->getElement($i,$j) * $v1->get($j);
  1153.             }
  1154.         }
  1155.         return new Math_Vector($data);
  1156.     }/*}}}*/
  1157.  
  1158.     // Static operations
  1159.  
  1160.     /**
  1161.      * Create a matrix from a file, using data stored in the given format
  1162.      *
  1163.      * Lines starting with '#' will be assumed to be comments and will be skipped
  1164.      *
  1165.      * @static
  1166.      * @access public
  1167.      * @param string $filename name of file containing matrix data
  1168.      * @param optional string $format one of 'serialized' (default) or 'csv'
  1169.      * @return object a Math_Matrix instance on success, a PEAR_Error otherwise
  1170.      */
  1171.     function &readFromFile ($filename, $format='serialized') {/*{{{*/
  1172.         if (!file_exists($filename) || !is_readable($filename)) {
  1173.             return PEAR::raiseError('File cannot be opened for reading');
  1174.         }
  1175.         if (filesize($filename) == 0) {
  1176.             return PEAR::raiseError('File is empty');
  1177.         }
  1178.         if ($format == 'serialized') {
  1179.             if (function_exists("file_get_contents")) {
  1180.                 $objser = file_get_contents($filename);
  1181.             } else {
  1182.                 $objser = implode("",file($filename));
  1183.             }
  1184.             $obj = unserialize($objser);
  1185.             if (Math_Matrix::isMatrix($obj)) {
  1186.                 return $obj;
  1187.             } else {
  1188.                 return PEAR::raiseError('File did not contain a Math_Matrix object');
  1189.             }
  1190.         } else { // assume CSV data
  1191.             $data = array();
  1192.             $lines = file($filename);
  1193.             foreach ($lines as $line) {
  1194.                 if (preg_match('/^#/', $line)) {
  1195.                     continue;
  1196.                 } else {
  1197.                     $data[] = explode(',',trim($line));
  1198.                 }
  1199.             }
  1200.             $m =& new Math_Matrix();
  1201.             $e = $m->setData($data);
  1202.             if (PEAR::isError($e)) {
  1203.                 return $e; 
  1204.             } else {
  1205.                 return $m;
  1206.             }
  1207.         }
  1208.     }/*}}}*/
  1209.  
  1210.     /**
  1211.      * Writes matrix object to a file using the given format
  1212.      *
  1213.      * @static
  1214.      * @access public
  1215.      * @param object Math_Matrix $matrix the matrix object to store
  1216.      * @param string $filename name of file to contain the matrix data
  1217.      * @param optional string $format one of 'serialized' (default) or 'csv'
  1218.      * @return mixed TRUE on success, a PEAR_Error otherwise
  1219.      */
  1220.     function writeToFile (&$matrix, $filename, $format='serialized') {/*{{{*/
  1221.         if (!Math_Matrix::isMatrix($matrix)) {
  1222.             return PEAR::raiseError("Parameter must be a Math_Matrix object");
  1223.         }
  1224.         if ($matrix->isEmpty()) {
  1225.             return PEAR::raiseError("Math_Matrix object is empty");
  1226.         }
  1227.         if ($format == 'serialized') {
  1228.             $data = serialize($matrix);
  1229.         } else {
  1230.             $data = '';
  1231.             list($nr, $nc) = $matrix->getSize();
  1232.             for ($i=0; $i<$nr; $i++) {
  1233.                 $row = $matrix->getRow($i);
  1234.                 if (PEAR::isError($row)) {
  1235.                     return $row;
  1236.                 }
  1237.                 $data .= implode(',', $row)."\n";
  1238.             }
  1239.         }
  1240.         $fp = fopen($filename, "w");
  1241.         if (!$fp) {
  1242.             return PEAR::raiseError("Cannot write matrix to file $filename");
  1243.         }
  1244.         fwrite($fp, $data);
  1245.         fclose($fp);
  1246.         return true;
  1247.     }/*}}}*/
  1248.  
  1249.     /**
  1250.      * Checks if the object is a Math_Matrix instance
  1251.      *
  1252.      * @static
  1253.      * @access public
  1254.      * @param object Math_Matrix $matrix
  1255.      * @return boolean TRUE on success, FALSE otherwise
  1256.      */
  1257.     function isMatrix (&$matrix) {/*{{{*/
  1258.         if (function_exists("is_a")) {
  1259.             return is_a($matrix, "Math_Matrix");
  1260.         } else {
  1261.             return (get_class($matrix) == "math_matrix");
  1262.         }
  1263.     }/*}}}*/
  1264.  
  1265.     /**
  1266.      * Returns a Math_Matrix object of size (nrows, ncols) filled with a value
  1267.      *
  1268.      * @static
  1269.      * @access public
  1270.      * @param integer $nrows number of rows in the generated matrix
  1271.      * @param integer $ncols number of columns in the generated matrix
  1272.      * @param numeric $value the fill value
  1273.      * @return object a Math_Matrix instance on success, PEAR_Error otherwise
  1274.      */
  1275.     function &makeMatrix ($nrows, $ncols, $value) {/*{{{*/
  1276.         if (!is_int($nrows) && is_int($ncols) && !is_numeric($value)) {
  1277.             return PEAR::raiseError('Number of rows, columns, and a numeric fill value expected');
  1278.         }
  1279.         for ($i=0; $i<$nrows; $i++) {
  1280.             $m[$i] = explode(":",substr(str_repeat($value.":",$ncols),0,-1));
  1281.         }
  1282.         return new Math_Matrix($m);
  1283.     }/*}}}*/
  1284.  
  1285.     /**
  1286.      * Returns the Math_Matrix object of size (nrows, ncols), filled with the value 1 (one)
  1287.      *
  1288.      * @static
  1289.      * @access public
  1290.      * @param integer $nrows number of rows in the generated matrix
  1291.      * @param integer $ncols number of columns in the generated matrix
  1292.      * @return object a Math_Matrix instance on success, PEAR_Error otherwise
  1293.      * @see Math_Matrix::makeMatrix()
  1294.      */
  1295.     function &makeOne ($nrows, $ncols) {/*{{{*/
  1296.         return Math_Matrix::makeMatrix ($nrows, $ncols, 1);
  1297.     }/*}}}*/
  1298.  
  1299.     /**
  1300.      * Returns the Math_Matrix object of size (nrows, ncols), filled with the value 0 (zero)
  1301.      *
  1302.      * @static
  1303.      * @access public
  1304.      * @param integer $nrows number of rows in the generated matrix
  1305.      * @param integer $ncols number of columns in the generated matrix
  1306.      * @return object a Math_Matrix instance on success, PEAR_Error otherwise
  1307.      * @see Math_Matrix::makeMatrix()
  1308.      */
  1309.     function &makeZero ($nrows, $ncols) {/*{{{*/
  1310.         return Math_Matrix::makeMatrix ($nrows, $ncols, 0);
  1311.     }/*}}}*/
  1312.  
  1313.     /**
  1314.      * Returns a square unit Math_Matrix object of the given size
  1315.      *
  1316.      * A unit matrix is one in which the elements follow the rules:
  1317.      *  e[i][j] = 1, if i == j
  1318.      *  e[i][j] = 0, if i != j
  1319.      * Such a matrix is also called an 'identity matrix'
  1320.      *
  1321.      * @static
  1322.      * @access public
  1323.      * @param integer $size number of rows and columns in the generated matrix
  1324.      * @return object a Math_Matrix instance
  1325.      * @see Math_Matrix::makeIdentity()
  1326.      */
  1327.     function &makeUnit ($size) {/*{{{*/
  1328.         for ($i=0; $i<$size; $i++) {
  1329.             for ($j=0; $j<$size; $j++) {
  1330.                 if ($i == $j) {
  1331.                     $data[$i][$j] = (float) 1.0;
  1332.                 } else {
  1333.                     $data[$i][$j] = (float) 0.0;
  1334.                 }
  1335.             }
  1336.         }
  1337.         return new Math_Matrix($data);
  1338.     }/*}}}*/
  1339.  
  1340.     /**
  1341.      * Returns the identity matrix of the given size. An alias of Math_Matrix::makeUnit()
  1342.      *
  1343.      * @static
  1344.      * @access public
  1345.      * @param integer $size number of rows and columns in the generated matrix
  1346.      * @return object a Math_Matrix instance
  1347.      * @see Math_Matrix::makeUnit()
  1348.      */
  1349.     function &makeIdentity($size) {/*{{{*/
  1350.         return Math_Matrix::makeUnit($size);
  1351.     }/*}}}*/
  1352.  
  1353.     /**
  1354.      * Solves a system of linear equations: Ax = b
  1355.      *
  1356.      * A system such as:
  1357.      * <pre>
  1358.      *     a11*x1 + a12*x2 + ... + a1n*xn = b1
  1359.      *     a21*x1 + a22*x2 + ... + a2n*xn = b2
  1360.      *     ...
  1361.      *     ak1*x1 + ak2*x2 + ... + akn*xn = bk
  1362.      * </pre>
  1363.      * can be rewritten as:
  1364.      * <pre>
  1365.      *     Ax = b
  1366.      * </pre>
  1367.      * where: 
  1368.      * - A is matrix of coefficients (aij, i=1..k, j=1..n), 
  1369.      * - b a vector of values (bi, i=1..k),
  1370.      * - x the vector of unkowns (xi, i=1..n)
  1371.      * Using: x = (Ainv)*b
  1372.      * where:
  1373.      * - Ainv is the inverse of A
  1374.      *
  1375.      * @static
  1376.      * @access public
  1377.      * @param object Math_Matrix $a the matrix of coefficients
  1378.      * @param object Math_Vector $b the vector of values
  1379.      * @return mixed a Math_Vector object on succcess, PEAR_Error otherwise
  1380.      * @see vectorMultiply()
  1381.      */
  1382.     function solve($a, $b) {
  1383.         // check that the vector classes are defined
  1384.         if (!Math_Matrix::isMatrix($a) && !Math_VectorOp::isVector($b)) {
  1385.             return PEAR::raiseError('Incorrect parameters, expecting a Math_Matrix and a Math_Vector');
  1386.         }
  1387.         $e = $a->invert();
  1388.         if (PEAR::isError($e)) {
  1389.             return $e;
  1390.         }
  1391.         return $a->vectorMultiply($b);
  1392.     }
  1393.  
  1394.     /**
  1395.      * Solves a system of linear equations: Ax = b, using an iterative error correction algorithm
  1396.      *
  1397.      * A system such as:
  1398.      * <pre>
  1399.      *     a11*x1 + a12*x2 + ... + a1n*xn = b1
  1400.      *     a21*x1 + a22*x2 + ... + a2n*xn = b2
  1401.      *     ...
  1402.      *     ak1*x1 + ak2*x2 + ... + akn*xn = bk
  1403.      * </pre>
  1404.      * can be rewritten as:
  1405.      * <pre>
  1406.      *     Ax = b
  1407.      * </pre>
  1408.      * where: 
  1409.      * - A is matrix of coefficients (aij, i=1..k, j=1..n), 
  1410.      * - b a vector of values (bi, i=1..k),
  1411.      * - x the vector of unkowns (xi, i=1..n)
  1412.      * Using: x = (Ainv)*b
  1413.      * where:
  1414.      * - Ainv is the inverse of A
  1415.      *
  1416.      * The error correction algorithm uses the approach that if:
  1417.      * - xp is the approximate solution
  1418.      * - bp the values obtained from pluging xp into the original equation
  1419.      * We obtain
  1420.      * <pre>
  1421.      *     A(x - xp) = (b - bp),
  1422.      *     or
  1423.      *     A*xadj = (b - bp)
  1424.      * </pr>
  1425.      * where:
  1426.      * - xadj is the adjusted value (= Ainv*(b - bp))
  1427.      * therefore, we calculate iteratively new values of x using the estimated
  1428.      * xadj and testing to check if we have decreased the error.
  1429.      *
  1430.      * @static
  1431.      * @access public
  1432.      * @param object Math_Matrix $a the matrix of coefficients
  1433.      * @param object Math_Vector $b the vector of values
  1434.      * @return mixed a Math_Vector object on succcess, PEAR_Error otherwise
  1435.      * @see vectorMultiply()
  1436.      * @see invert()
  1437.      * @see Math_VectorOp::add()
  1438.      * @see Math_VectorOp::substract()
  1439.      * @see Math_VectorOp::length()
  1440.      */
  1441.     function solveEC($a, $b) {/*{{{*/
  1442.         $ainv = $a->clone();
  1443.         $e = $ainv->invert();
  1444.         if (PEAR::isError($e)) {
  1445.             return $e;
  1446.         }
  1447.         $x = $ainv->vectorMultiply($b);
  1448.         if (PEAR::isError($x)) {
  1449.             return $x;
  1450.         }
  1451.         // initial guesses
  1452.         $bprime = $a->vectorMultiply($x);
  1453.         if (PEAR::isError($bprime)) {
  1454.             return $bprime;
  1455.         }
  1456.         $err = Math_VectorOp::substract($b, $bprime);
  1457.         $adj = $ainv->vectorMultiply($err);
  1458.         if (PEAR::isError($adj)) {
  1459.             return $adj;
  1460.         }
  1461.         $adjnorm = $adj->length();
  1462.         $xnew = $x;
  1463.         
  1464.         // compute new solutions and test for accuracy
  1465.         // iterate no more than 10 times
  1466.         for ($i=0; $i<10; $i++) {
  1467.             $xnew = Math_VectorOp::add($x, $adj);
  1468.             $bprime = $a->vectorMultiply($xnew);
  1469.             $err = Math_VectorOp::substract($b, $bprime);
  1470.             $newadj = $ainv->vectorMultiply($err);
  1471.             $newadjnorm = $newadj->length();
  1472.             // did we improve the accuracy?
  1473.             if ($newadjnorm < $adjnorm) {
  1474.                 $adjnorm = $newadjnorm;
  1475.                 $x = $xnew;
  1476.                 $adj = $newadj;
  1477.             } else { // we did improve the accuracy, so break;
  1478.                 break;
  1479.             }
  1480.         }
  1481.         return $x;
  1482.     }/*}}}*/
  1483.    
  1484. } // end of Math_Matrix class /*}}}*/
  1485.  
  1486. // vim: ts=4:sw=4:et:
  1487. // vim6: fdl=1:
  1488.  
  1489. ?>
  1490.